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ABSTRACT 


The  objective  of  this  research  was  to  develop  an  optimization  technique  that  can  be 
used  interactively  by  design  engineers  to  approach  an  optimal  design  with  minimal 
computational  effort.  The  technique  can  be  applied  to  both  continuous  and  discrete  values 
of  design  variables.  A  large  number  of  design  variables  can  be  also  considered. 

In  order  to  meet  the  objective,  an  optimization  procedure  was  developed  by  coupling 
the  finite  element  analysis  (FEA)  to  the  orthogonal  array  experimentation  technique,  because 
FEA  is  a  common  analysis  tool  for  design  engineers.  From  the  results  of  the  FEA  and  an 
orthogonal  array,  an  average  Jacobian  matrix  was  constructed  that  showed  the  average 
overall  sensitivity  of  the  design  variables.  These  sensitivities  were  then  used  to  optimize  the 
design  parameters.  The  process  could  then  be  repeated  at  the  discretion  of  the  engineer  until 
a  satisfactory  design  is  obtained.  In  general,  the  designer  can  predict  and  control  the 
number  of  FEA  calculations  before  an  optimization  process  so  that  one  can  plan  a  budget 
and  time  for  an  optimal  design. 

Some  examples  of  structural  optimization  with  truss  and  frame  structures  with 
continuous  and  discrete  values  of  design  variables  were  studied  using  the  technique 
developed  in  this  paper.  Their  optimal  solutions  were  found  with  small  numbers  of 
iterations. 
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I.  INTRODUCTION 


Optimization  is  a  process  that  can  be  seen  in  almost  every  aspect  of  life.  Humans  are 
always  striving  to  go  faster,  climb  higher,  run  farther,  and  engineering  is  no  different. 
Engineering  requires  optimization  so  that  we  can  design  faster  planes  and  cars,  that  use  less 
fuel,  are  lighter,  stronger,  and  more  comfortable.  Engineers  have  been  optimizing  designs 
since  the  beginning,  but  recent  advances  in  computing  have  made  numerical  optimization 
techniques  a  more  effective  way  than  the  original  trial-and-error  and  experience  based 
optimization.  The  present  methods  based  on  numerical  search  techniques  are  effective  and 
easy  to  understand  and  implement  on  design  problems  that  don't  have  many  design 
variables  (Vanderplaats,  1984).  However,  it  is  difficult  to  estimate  the  computational  time 
and  budget  to  converge  to  an  optimal  design  using  most  of  the  numerical  optimization 
techniques  available  today  unless  the  design  engineer  has  experience  with  the  same  kind  of 
design  problem  in  advance.  The  computational  costs  increase  highly  nonlinearly  as  the 
number  of  design  variables  increases.  Therefore,  application  of  a  numerical  optimization 
technique  to  a  practical  engineering  design  has  been  hampered. 

The  objective  of  this  research  was  to  develop  an  efficient  method  for  achieving  the 
optimal  design  with  minimal  computational  effort.  Also,  the  developed  technique  allows  the 
engineer  to  interact  with  the  procedure  and  to  estimate  how  much  time  is  to  be  spent  on  the 
optimization  depending  on  the  number  of  design  variables  (it  does  not  matter  how  many 
design  variables  there  are),  and  how  close  to  the  actual  optimum  the  design  must  be, 
balancing  schedule  and  cost.  The  approach  was  to  develop  an  optimization  procedure  by 
combining  the  finite  element  method  and  orthogonal  arrays.  The  emphasis  was  in  structural 
mechanics,  but  could  be  applied  to  other  disciplines.  Because  finite  element  analyses  are 
time  consuming  in  general,  the  approach  was  to  minimize  the  number  of  iterations 
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associated  with  the  FEA.  This  led  to  the  incorporation  of  orthogonal  arrays  in  the 
methodology  in  order  to  reduce  the  number  of  FEA  computations. 

Orthogonal  arrays  were  developed  initially  for  discrete  values  of  design  variables. 
Therefore,  they  work  well  for  the  discrete  case.  Further,  the  optimization  technique  was 
also  applied  to  continuous  design  variables  with  design  constraints. 

The  goal  is  not  necessarily  to  get  the  actual  optimum  value,  assuming  that  one  exists, 
but  to  develop  a  more  practical  approach  that  will  achieve  a  result  close  to  the  optimum  with 
reasonable  effort.  In  addition  to  developing  a  methodology  that  would  lead  to  an  optimal 
design,  it  was  necessary  to  identify  the  overall  sensitivity  of  the  various  parameters.  This 
information  provides  an  engineer  with  further  physical  understanding  of  the  design 
problem. 

Although  the  technique  may  not  achieve  the  most  optimal  value,  it  will  give  the 
designer  the  value  of  parameters  close  to  the  optimum  with  reasonable  effort  and  drastically 
reduced  computational  time.  Just  like  if  you  want  a  higher  quality  car  you  might  buy  a 
luxury  car,  or  a  faster  car  you  might  buy  a  sports  car,  you  must  pay  the  price.  Similarly,  if 
you  want  the  most  optimum  design  you  must  pay  for  it  with  either  higher  design  and 
manufacturing  costs  or  more  computer  time.  On  the  other  hand,  if  you  are  willing  to  accept 
a  good  solid  car  you  might  buy  a  basic  family  car.  This  is  similar  to  spending  a  reasonable 
amount  of  time  optimizing  and  getting  a  close  to  optimal  design  without  the  high  computer 
computational  times.  The  latter  approach  is  therefore  the  goal  which  should  lead  to  a  more 
practical  optimization  approach. 

This  paper  will  first  discuss  the  background  of  the  Taguchi  Method  and  the 
orthogonal  arrays  that  he  has  incorporated  into  his  design  of  experiments  methodology. 
This  will  be  followed  by  a  brief  technical  description  of  the  orthogonal  arrays  as  described 
by  Peace  (1993).  Next  will  be  the  method  and  procedure  developed  in  this  research.  The 
final  section  will  include  some  examples  of  the  technique. 
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II.  TAGUCHI  METHOD  AND  ORTHOGONAL  ARRAYS 


The  optimization  approach  developed  in  this  research  centered  around  the  concepts 
developed  by  Taguchi  for  quality  engineering  and  then  coupling  his  concepts  with  the  finite 
element  method  to  achieve  an  optimal  design  methodology.  The  following  brief  description 
of  orthogonal  arrays  can  be  found  in  Peace  (Peace,  1993). 

The  orthogonal  array  was  the  foundation  of  the  methodology  that  Taguchi  developed 
for  designing  an  experiment.  Although  other  methods  of  conducting  an  experiment  are 
possible,  like  a  full  factorial  approach,  the  orthogonal  array  has  traditionally  been  associated 
with  the  Taguchi  experimentation  technique.  The  orthogonal  array  is  repeatable,  cost 
effective,  and  efficient  in  that  only  a  small  amount  of  data  must  be  collected  and  then  that 
data  can  easily  be  translated  into  meaningful  and  verifiable  conclusions  (Peace,  1993). 

The  orthogonal  array  was  originally  developed  by  Fisher  of  England.  It  was 
developed  during  his  efforts  to  control  error  in  an  experiment  (Ealey,  1988).  Later,  Taguchi 
adapted  the  orthogonal  array  in  order  to  measure  the  effect  of  a  factor  being  studied  on  the 
average  result  (Peace,  1993). 

The  details  of  the  orthogonal  array  will  follow,  but  the  concept  of  orthogonality 
refers  to  the  statistically  independent  or  balanced  parameters  that  make  up  the  columns  of  the 
orthogonal  array.  One  major  objective  of  this  research  was  the  application  of  the  technique 
in  both  discrete  and  continuous  parameter  problems.  However,  the  orthogonal  array 
approach  lends  itself  directly  only  to  problems  with  discrete  parameters.  For  problems 
associated  with  discrete  parameters,  the  orthogonal  array  is  used  directly.  But,  for  problems 
with  continuous  variables,  this  research  will  show  that  taking  the  minimum  and  maximum 
values  in  the  range  of  each  of  the  parameters,  and  then  putting  those  values  into  the 
orthogonal  zirray,  is  an  effective  way  of  dealing  with  continuous  variables.  A  discussion  of 
how  this  was  dealt  with  will  follow  this  section. 
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In  a  typical  orthogonal  array  as  shown  in  Table  1  (Taguchi,  1987),  each  level  (here, 
1  or  2)  has  an  equal  number  of  occurrences  in  each  column.  Also  notice  the  relationship 
between  one  column  and  another.  For  each  level  within  one  column,  each  level  within  any 
other  column  will  occur  an  equal  number  of  times  as  well.  This  introduces  the  statistical 
independence  or  balance  into  the  orthogonal  array.  Since  each  column  is  orthogonal  to  the 
others,  if  the  results  associated  with  one  level  of  a  specific  factor  are  much  different  at 
another  level,  it  is  because  changing  that  factor  from  one  level  to  another  has  a  strong  impact 
on  the  design  parameter  being  measured.  Because  the  levels  of  the  other  factors  are 
occurring  an  equal  number  of  times,  any  effect  by  these  other  factors  should  be  canceled  out 
(Peace,  1993). 


Table  1.  L8(2'^)  Orthogonal  Array 


Parameter  No. 

No.  of  Runs 

1 

2 

3 

4 

5 

6 

7 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

2 

2 

2 

2 

3 

1 

2 

2 

1 

1 

2 

2 

4 

1 

2 

2 

2 

2 

1 

1 

5 

2 

1 

2 

1 

2 

1 

2 

6 

2 

1 

2 

2 

1 

2 

1 

7 

2 

2 

1 

1 

2 

2 

1 

8 

2 

2 

1 

2 

1 

1 

2 

The  design  of  an  orthogonal  array  experiment,  although  balanced,  does  not  require 
that  all  combinations  of  all  factors  be  tested.  The  full  factorial  experiment  studying  7  factors 
would  require  128  (=  2”^)  experimental  runs,  compared  with  the  8  runs  in  the  orthogonal 
array.  The  advantage  of  the  orthogonal  array  can  be  seen  in  its  simplicity,  time  savings,  and 
therefore  cost  savings. 
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It  should  be  noted,  that  the  columns  in  the  orthogonal  array  can  represent  not  only 
design  factors,  but  also  interactions.  Therefore,  by  using  sound  engineering  judgment,  the 
orthogonal  array  can  include  critical  interactions.  The  orthogonal  array,  does  not  study  all 
interaction  combinations  like  a  full  factorial  experiment,  but  the  balanced  nature  of  the 
orthogonal  array  can  compensate  for  this. 

The  terminology  of  the  orthogonal  array  La(BC)  is  as  follows.  The  subscript  of  L, 
which  here  is  A,  represents  the  number  of  experimental  runs  which  can  be  conducted  in  the 
experiment.  B  denotes  the  number  of  levels  within  each  column,  and  C  denotes  the  number 
of  columns  (factors  and/or  interactions  that  can  be  studied)  in  the  orthogonal  array. 
Therefore,  the  above  orthogonal  array  LgCl^)  contains  eight  experimental  runs.  Each 
column  contains  '2'  levels,  meaning  each  factor  will  posses  two  levels.  The  7'  means  that 
up  to  seven  factors  and  interactions  can  be  incorporated  into  the  experiment. 
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III.  METHOD  &  PROCEDURE 


As  discussed  earlier,  the  procedure  that  has  been  developed  focuses  on  the 
orthogonal  array.  The  following  example  is  written  for  a  structural  optimization  problem 
with  stress,  displacement,  nahiral  frequencies,  etc.  as  the  design  constraints,  but  in  general 
could  be  applied  to  other  disciplines.  The  design  variables  varied  from  a  minimum  value  to 
a  maximum  value  continuously  or  discretely.  The  objective  was  to  reduce  the  weight  of  the 
structure  supporting  a  load  while  not  exceeding  the  design  constraints  in  any  of  the  members 
of  the  structure.  In  this  example,  the  cross-sectional  areas  are  considered  as  the  design 
variables,  the  weight  of  the  structure  was  the  objective  function,  and  the  stresses, 
displacements,  and  natural  frequencies  were  the  constraints. 

The  general  procedure  of  the  optimization  technique  that  was  developed  begins  by 
constructing  the  proper  orthogonal  array  for  the  given  design  variables  using  the  entire  range 
of  values  for  each  variable.  If  the  design  variables  are  continuous,  then  the  values  of  1  and 
2  are  used  in  the  orthogonal  array  corresponding  to  the  minimum  and  maximum  values 
respectively.  If  the  design  variables  are  discrete,  then  the  lowest  discrete  value  is  used  for 
level  1  in  the  orthogonal  array,  the  next  lowest  value  is  used  for  level  2  in  the  orthogonal 
array,  etc. 

The  next  step  is  to  use  the  design  values  in  each  row  of  the  orthogonal  array  and  run 
a  FEM  program  (Kwon,  1997)  for  each  row  using  the  corresponding  design  values  in  the 
columns  of  the  orthogonal  array.  The  inputs  to  the  FEM  code  are  the  cross-sectional  areas, 
and  the  output  from  the  FEM  code  are  the  values  of  normalized  stress,  displacements, 
natural  frequencies,  etc.,  P-.  Here  P-  is  the  constraint  function  and  is  defined  as 

p  M<i  (1) 

'  5,. 
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where  /,  is  a  computed  value  from  the  finite  element  analysis  such  as  stress,  displacement, 
or  natural  frequency.  And,  5,  is  the  allowable  value  of  /,  such  as  failure  strength, 
displacement,  or  natural  frequency,  and  is  a  factor  of  safety  ( Jfc,.  >  1 ) . 

After  all  cases  have  been  run  in  the  orthogonal  array,  the  results  are  plotted  vs. 

which  is  the  sensitivity  of  a  normalized  constraint  value  with  respect  to  a  design  variable, 

AP 

namely,  the  average  overall  sensitivities  (slopes)  — -  are  calculated.  With  the  desired  value 

Axj 

for  all  P.  =  1,  a  target  value  of  P.  was  selected  using  a  one-dimensional  (1-D)  search 

technique  in  the  direction  determined  by  the  average  overall  sensitivities  until  a  constraint 
boundary  was  encountered. 

The  next  step  is  to  take  the  best  value  from  all  of  the  orthogonal  array  runs  that  meet 
the  design  criteria  and  constraints,  although  it  will  be  shown  that  the  final  solution  is 
relatively  independent  of  the  starting  point.  It  is  recommended  that  the  starting  point  be  in 
the  feasible  domain  for  all  parameters,  because  the  technique  is  attempting  to  approach  the 
optimum  design  from  within  the  feasible  domain.  If  no  cases  in  the  orthogonal  array  results 
are  feasible,  the  most  conservative  design  may  be  used  as  the  initial  design.  These  cross- 
sectional  areas  will  then  form  the  vector,  {xi}.  These  values  are  then  used  to  run  the  FEA 
code,  which  produces  a  vector  of  P- values  called  {Pi}.  The  superscript  denotes  the  iteration 
number. 

From  the  orthogonal  array,  the  average  overall  sensitivity  of  each  parameter  over  its 
entire  domain  is  calculated.  This  is  slightly  different  from  the  common  optimization 
approach  of  calculating  the  Jacobian  (which  contains  the  local  sensitivities  calculated  at  one 
point).  In  the  traditional  approach,  the  Jacobian  matrix  is  determined  with  a  fixed  value  of 
each  design  variable.  The  Jacobian  matrix  (Lang,  1968)  is 
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1 _ 

dP, 

dP, 

dx^ 

dx„ 

dP^ 

dP2 

dx^ 

dx2 

dx„ 

dK 

dK 

dxi 

dx2 

dx 

(2) 


where  n  is  the  number  of  design  variables  and  m  is  the  number  of  constraints.  In  the 
present  study,  the  averaged  Jacobian  matrix  over  the  domain  of  each  variable  can  be 
assembled  from  the  average  overall  sensitivities  instead  of  being  computed  at  a  specific  set 
of  design  values.  Thus  the  averaged  Jacobian  matrix  is 


[^1 


Ax. 

Ax, 

Ax 

AP, 

Ax, 

Ax, 

Ax„ 

Ax, 

Ax, 

Ax„  ^ 

(3) 


In  general ,  the  off-diagonal  terms  in  the  first  n  rows  are  small  compared  with  the 
diagonal  terms  and  are  therefore  neglected.  Thus  [7]  can  be  rewritten  as 


9 


0 


■  ^ 
Ax, 

0 


0 

Ax2 


0 


0 


AP. 


n+l 


Ax, 


0 

AX2 


AP. 


n+l 


Ax. 


AP„  AP„  AP„ 

Ax,  Ax2  Ax„ 


(4) 


Out  of  the  terms  in  the  averaged  Jacobian  matrix,  the  absolute  maximum  in  each  column  is 
used  because  that  value  is  corresponds  to  the  most  constraining  value  when  the  design 
variables  are  updated.  Using  these  sensitivities,  the  following  equation  is  solved  to  compute 
the  change  of  the  design  variables 


AP,= 


K^jJ 


Axj 


(5) 


for  i=l...n  and  j=l...m.  Here  Ax,  is  simply  the  needed  change  in  each  design  variable  to 
achieve  an  optimal  value  which  is  to  be  solved  for  using 


Tap.^  * 


where 


(7) 
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Here,  AP.  is  the  distance  to  be  traveled  in  the  direction  determined  by  the  orthogonal  array 
from  the  current  value  of  P.  to  the  target  value  P..  The  target  value  P.  is  the  value  to  be 

determined  from  the  1-D  search.  This  value  can  be  found  by  extrapolation  until  a  constraint 
is  violated  and  then  interpolated  until  the  constraint  boundary  is  found.  Many  1-D  search 
techniques  exist  that  can  efficiently  accomplish  this  task  such  as  the  bisection  method  or  the 
golden  section  method  (Vanderplaats,  1984). 

Axj  is  then  added  to  the  previous  design  value  xj^  to  yield  the  new  design  values  xj^. 
The  new  design  values  are  then  used  in  running  the  FEA  program  again  to  verify  the  result. 

At  this  point,  if  any  or  all  of  the  values  in  Pj^  are  close  to  1,  then  this  is  assumed  to 
be  very  close  to  the  optimum  design,  and  the  analysis  can  be  stopped,  with  the  design 
values  becoming  xj^.  If  further  reduction  is  required,  then  the  range  of  the  parameter  values 
is  refined  by  using  ±5%  of  the  values  in  Xj2.  The  procedure  is  then  repeated  using  the  new 
refined  maximum  and  minimum  values  as  the  limits  in  the  orthogonal  array  until  a 
satisfactory  design  is  obtained.  Usually,  more  than  one  set  of  orthogonal  arrays  is  not 
necessary  unless  a  very  refined  optimal  design  is  required.  This  further  refinement  is  the 
topic  of  follow-on  studies. 

For  a  discrete  value  of  a  design  variable,  the  change  of  design  variable  Axj  is 
compared  to  the  allowable  change  of  the  design  variable.  If  the  magnitude  of  the  calculated 
Axj  is  greater  than  the  magnitude  of  the  allowable  change,  then  Axj  is  set  equal  to  the 
allowable  value.  Otherwise,  the  design  variable  does  not  change.  This  process  continues 
until  there  is  no  longer  any  change  in  the  design. 

Figure  1  outlines  the  procedure  described  above. 
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Optimization  Procedure  using  Orthogonal  Arrays _ 

1 .  Construct  the  Proper  Orthogonal  Array. 

2 .  Select  an  Initial  (Feasible)  Design. 

3 .  Calculate  the  Average  Overall  Sensitivities. 

4 .  Update  the  Design  Variables  using  a  1-D  Search  Technique. 

5 .  If  the  design  is  satisfactory,  then  stop  the  procedure.  Otherwise,  repeat 

the  procedure  with  a  refined  domain  for  further  optimization. _ 

Figure  1.  Procedure  Summary 


12 


IV.  RESULTS  &  DISCUSSION 


In  this  section,  three  examples  will  be  shown  to  illustrate  the  technique  described  in 
the  previous  section.  The  first  example  is  a  statically  independent  3-bar  truss  with 
continuous  design  variables  and  stress/displacement  constraints.  The  second  example  is  a 
statically  independent  6-bar  truss  in  a  star  shape  with  discrete  design  variables  and  only 
stress  constraints.  The  third  example  is  a  12-bar  frame  structure  with  continuous  design 
variables.  It  has  stress,  maximum  displacement,  and  minimum  natural  frequency 
constraints.  The  frame  problem  in  the  third  example  is  first  solved  using  a  conventional 
numerical  optimization  technique.  Sequential  Quadratic  Programming  (SQP).  Next,  the 
frame  problem  is  solved  with  orthogonal  arrays  using  the  techniques  developed  in  this 
research.  Finally,  the  results  of  the  two  optimization  techniques  are  compared. 

A.  STATICALLY  INDETERMINATE  THREE-BAR  TRUSS  EXAMPLE 


This  first  example  is  a  statically  independent  three-bar  truss  structure  subjected  to  a 
load  of  50,000  lbs  to  the  right  and  100,000  lbs  down.  The  objective  is  to  minimize  the 
weight  of  the  structure  such  that  each  member  does  not  exceed  the  failure  strength  of  1 18  ksi 
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and  the  displacement  at  the  node  where  the  load  is  applied  should  not  exceed  0.12  inches 
horizontally  and  -0.06  inches  vertically.  A  factor  of  safety  of  1.5  was  used  for  the  stress 
and  displacement  values.  Initially,  the  minimum  cross-sectional  area  is  0.25  in^,  and  the 
maximum  cross-sectional  area  is  2.00  in^  for  all  members.  The  design  variables  (cross- 
sectional  area)  in  this  problem  were  continuous.  These  areas  are  then  put  into  the  following 
L4(23)  orthogonal  array  (Phadke,  1989)  (See  Table  2). 


Table  2.  14(2^)  Orthogonal  Array 


Run# 

Area  1 

Area  2 

Area  3 

1 

0.25 

0.25 

2 

0.25 

2.00 

2.00 

3 

2.00 

0.25 

2.00 

4 

2.00 

2.00 

0.25 

Table  3  shows  the  results  from  running  the  FEM  code  for  each  of  the  4  rows  in  the 
orthogonal  array,  Table  2. 


Table  3.  FEA  Results 


Run# 

PI 

P2 

P3 

P4 

P5 

wf  1 

1 

1.03 

1.69 

4.33 

1.53 

6.43 

2 

0.27 

0.32 

0.49 

0.41 

0.73 

3 

0.28 

0.35 

0.62 

0.41 

0.91 

78.71 

4 

-0.15 

0.52 

3.22 

-0.22 

4.77 

104.67 

The  first  three  columns  (PI  -  P3)  are  the  normalized  stresses  and  the  next  two 
columns  (P4  -  P5)  are  the  normalized  displacements  of  the  node  where  the  load  is  applied. 

Figure  3  is  a  sensitivity  plot  for  the  normalized  stress  in  element  1  with  respect  to  the 
three  different  cross-sectional  areas  (Al,  A2,  A3).  Similarly,  Figure  4  is  a  sensitivity  plot 
for  the  normalized  stress  in  element  2  with  respect  to  the  three  different  cross-sectional 
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areas.  And  finally,  Figure  5  is  a  sensitivity  plot  for  the  normalized  stress  in  element  3  with 
respect  to  the  three  different  cross-sectional  areas.  In  all  these  three  plots,  the  dotted  lines 
represent  the  stress  relative  to  element  1,  the  dashed  lines  represent  the  stress  relative  to 
element  2,  and  the  dash-dot  lines  represent  the  stress  relative  to  element  3.  From  these 
plots,  one  can  see  that  the  slope  of  the  normalized  stress  in  each  member  vs.  the  area  of  that 
same  member  is  the  steepest  (or  close  to  the  steepest)  and  is  therefore  has  the  highest 
average  overall  sensitivity. 

We  then  take  only  the  sensitivity  plot  for  element  1  with  respect  to  the  cross- 
sectional  area  of  element  1  and  plot  this  sensitivity  in  Figure  6.  Similarly,  Figure  7  is  the 
sensitivity  of  element  2  with  respect  to  the  cross-sectional  area  in  element  2.  Finally,  Figure 
8  is  for  element  3.  Notice  that  the  slope  in  Figure  8  is  the  steeper  than  Figures  6  and  7, 
indicating  that  the  stress  level  in  element  3  is  much  more  sensitive  to  changes  in  cross- 
sectional  area  3  than  in  elements  1  and  2. 

Figures  9-11  show  the  sensitivities  of  the  normalized  horizontal  displacement,  P4, 
to  the  three  cross-sectional  areas.  Here  we  see  that  all  of  the  slopes  are  about  the  same 
indicating  that  the  horizontal  displacement  is  affected  equally  by  changes  in  all  three 
elements. 

Figures  12-14  show  the  sensitivities  of  the  normalized  horizontal  displacement,  P4, 
to  the  three  cross-sectional  areas.  Figure  14  shows  the  steepest  slope,  indicating  that 
changes  in  the  cross-sectional  area  of  element  3  will  have  the  greatest  impact  on  the  vertical 
displacement. 
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Umax 


In  this  analysis  we  will  use  run  3  as  the  starting  point  because  it  has  the  lowest 
weight  of  the  runs  where  no  constraints  are  violated.  The  starting  point  { }  is 

{a'}''  =  {2.00  0.25  2.00}  (8) 

with  the  resulting 

{P‘}^  =  {0.28  0.35  0.62  0.41  0.91}  (9) 

It  meets  the  criteria  that  none  of  the  constraints  are  violated,  therefore  the  design  variables 
are  all  in  the  feasible  domain.  This  run  results  in  a  weight  of  78.71  lbs.  This  becomes  the 
{A^ } vector  with  the  corresponding  result,  {pi}. 

In  order  to  optimize  the  design  variables,  the  following  calculations  were  performed 
for  each  member.  For  example,  for  member  1  equation  (5)  is  applied  for  both  stress  and 
displacement  constraints  with  a  target  value  of  0.9.  From  this,  three  values  (1  from  stress 
constraint  and  2  from  the  displacement  constraints)  of  AA  were  calculated,  then  the  largest 
absolute  value  was  chosen  to  update  the  design  variables.  The  resulting  design  variables 
were 

{a2}’'={0.75  0.25  1.90}  (10) 

In  order  to  verify  that  all  of  the  design  constraints  were  met,  the  FEM  code  was  run  using 
these  design  variables,  with  the  following  results 

{P^}^={0.65  0.65  0.63  0.97  0.94}  (11) 
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The  corresponding  weight  for  this  design  is  47.55  lbs,  a  considerable  weight  reduction  after 
only  five  FEA  computations.  The  results  show  that  the  new  design  variables  satisfy  the 
constraints  and  are  therefore  an  acceptable  solution.  However,  if  the  engineer  wants  to 
further  optimize  the  design  variables,  then  the  procedure  can  be  repeated  with  a  refined 
domain. 

The  refining  of  the  domain  is  done  individually  for  each  cross-sectional  area. 
{Amin}  is  calculated  from  the  equation 


similarly,  {  Amax)  is  calculated  from  the  equation 

K„}  =  {A'-'}+^{(l-|P'-i)K-'}}  (13) 

where  P  is  a  scaling  factor  to  adjust  the  refined  domain.  Once  again,  the  upper  and  lower 
values  are  checked  to  insure  that  the  values  remain  in  the  feasible  domain. 

In  general,  the  procedure  using  the  orthogonal  array  only  needs  to  be  carried  out 
once  in  most  practical  applications.  This  results  in  a  reasonably  optimum  design  without 
exhaustive  computations.  The  final  results  of  this  procedure  are  relatively  independent  of 
the  starting  point  in  the  analysis. 
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B.  STATICALLY  INDETERMINATE  STAR-SHAPED  TRUSS  EXAMPLE 


Figure  15.  6-Bar  Star-Shaped  Truss 


Following  the  procedure  outlined  in  the  previous  section,  this  example  shows  a 
statically  indeterminate  six-bar  tmss  structure  with  discrete  design  variables.  The  geometry 
was  chosen  to  insure  that  the  problem  was  fully  statically  indeterminate.  The  objective  here 
was  to  minimize  the  weight  of  the  structure,  while  not  exceeding  the  failure  stress  in  any  of 
the  six  members.  The  length  of  all  6  members  was  10  in.  The  discrete  cross-sectional  areas 
for  all  six  members  were  0.25  in^  ,  0.50  in^,  1.00  in^,  2.00  in^,  and  2.50  in^.  These 
represent  the  five  discrete  design  variables  to  be  used  in  the  orthogonal  array.  The  L25(5^) 
orthogonal  array  (Taguchi,  1987)  was  used  for  this  problem.  Table  4  shows  the  FEA 
results  for  each  run. 
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Table  4.  FEA  Results 


Run# 

PI 

P2 

P3 

P4 

P5 

P6 

WT 

1 

4.72 

-2.60 

3.26 

4.61 

3.57 

-1.04 

15.00 

2 

2.79 

-1.60 

1.85 

2.68 

2.15 

-0.53 

27.50 

3 

1.53 

-0.89 

1.00 

1.46 

1.20 

-0.27 

52.50 

4 

0.81 

-0.48 

0.52 

0.77 

0.63 

-0.14 

102.50 

5 

0.65 

-0.38 

0.42 

0.62 

0.51 

-0.11 

127.50 

6 

1.31 

-0.72 

0.90 

1.27 

0.99 

-0.29 

67.50 

7 

0.87 

-0.42 

0.68 

0.90 

0.62 

-0.28 

67.50 

8 

1.18 

-1.05 

0.33 

0.88 

1.17 

0.30 

67.50 

1.42 

-0.89 

0.85 

1.31 

1.15 

-0.17 

67.50 

1.65 

-0.63 

1.48 

1.80 

1.05 

-0.75 

67.50 

1.11 

-0.92 

0.40 

0.87 

1.06 

0.18 

72.50 

1.41 

-0.97 

0.75 

1.25 

1.20 

-0.05 

72.50 

0.96 

-0.52 

0.67 

0.94 

0.72 

-0.22 

72.50 

14 

1.12 

-0.14 

1.36 

1.43 

0.51 

-0.92 

72.50 

15 

0.96 

-0.48 

0.72 

0.97 

0.69 

-0.28 

72.50 

16 

0.79 

-0.40 

0.58 

0.79 

0.57 

-0.22 

82.50 

17 

1.04 

-0.83 

0.41 

0.84 

0.96 

0.13 

82.50 

18 

0.85 

-0.53 

0.52 

0.79 

0.68 

-0.11 

82.50 

19 

0.67 

-0.27 

0.59 

0.73 

0.44 

-0.29 

82.50 

20 

0.95 

-0.18 

1.07 

1.17 

0.47 

-0.69 

82.50 

21 

0.72 

-0.59 

0.26 

0.57 

0.68 

0.11 

87.50 

22 

0.61 

-0.27 

0.51 

0.65 

0.42 

-0.23 

87.50 

23 

0.85 

-0.26 

0.84 

0.98 

0.50 

-0.48 

87.50 

24 

0.91 

-0.48 

0.66 

0.91 

0.67 

-0.24 

87.50 

25 

0.72 

-0.35 

0.55 

0.73 

0.51 

-0.22 

87.50 

From  the  results  shown  in  Table  4,  the  overall  sensitivities  of  the  constraints  to  the 
AP 

design  variables  — -  between  each  level,  are  computed.  For  each  level  apply  Equation  5 
AA. 

and  then  compute  AAj .  If  AAj  is  less  than  the  amount  available  for  the  reduction  of  each 
design  variable,  then  the  value  of  the  design  variable  can  be  further  reduced. 

In  this  example,  the  7th  run  is  used  as  the  initial  design  and  is  updated  based  upon 
the  previously  described  procedure.  Notice  that  the  weight  of  the  7th  run  is  67.5  lbs., 
which  ends  up  being  very  close  to  the  final  value.  One  could  stop  here  and  get  a  very  good 
feasible  design.  The  final  result  was 
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{a2}’^  =  {0.50  0.25  1.00  2.00  2.00  0.25} 


(14) 


[P^Y  =  {0.91  -0.58  0.62  0.92  0.76  -0.16} 


(15) 


Because  each  of  the  values  in  the  results  vector  does  not  violate  any  of  the  design 
constraints,  the  values  in  Equation  14  are  acceptable  design  variables.  The  weight  for  this 
structure  is  60  lbs.,  a  reduction  of  40%  over  the  most  conservative  design  for  this  structure 
of  150  lbs  in  26  FEM  calculations. 


C.  OPTIMIZATION  OF  A  12-BAR  FRAME  STRUCTURE 

As  discussed  earlier,  the  12-bar  frame  structure  problem  is  first  solved  using 
common  numerical  optimization  techniques.  Next,  the  problem  is  solved  using  the 
procedures  developed  in  this  research.  Figure  16  shows  the  12-bar  frame  structure. 

1.  Solution  Using  Common  Numerical  Optimization  Techniques 

The  purpose  of  this  example  was  to  minimize  the  volume  of  a  frame  structure 
subject  to  a  horizontal  load  of  1000  lbs  at  the  top  of  the  left  edge,  and  1000  lbs.  at  the 
center  of  the  left  edge  (See  Figure  16).  The  constraints  on  the  problem  included  the 
maximum  stress  in  any  member  could  not  exceed  the  yield  strength  of  the  material,  the 
maximum  displacement  of  any  node  was  not  to  exceed  .125",  and  the  lowest  natural 
frequency  of  the  system  must  be  larger  than  7  Hz. 
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Figure  16.  12-Element  Frame  Structure 


The  cross-sectional  area  of  each  member  was  rectangular  with  a  2: 1  height-to-width 
ratio.  Therefore,  the  heights  of  each  of  the  members  made  up  the  12  design  variables.  The 
upper  limit  for  each  member  was  h=4.0",  and  the  lower  limit  for  each  member  was  h=0.1". 

The  main  program  executed  the  MATLAB  constr  optimization  routine  that  called  a 
'function  M-file'  (fun.m)  which  calculated  the  values  of  the  objective  and  constraints 
(Grace,  1992).  The  'function  M-file'  also  called  2  FEM  programs  that  were  used  to 
calculate  the  constraints.  The  first  FEM  program  {disp_frame_fem.m),  calculated  the  stress 
in  each  member  and  determined  the  maximum  horizontal  displacement  of  any  node.  The 
second  FEM  program  {eigjrame_fem.m),  calculated  the  natural  frequencies  of  the  system. 

The  problem  statement  for  the  constrained  optimization  problem  is  given  by  the 
following  objective  function,  F(h),  and  constraints,  g(h): 

12 

minimize:  F(h)-^bjh-L  (16) 

i=l 
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subject  to:  g(i)  =  - 1  <  0  for  i  = 

^yieU 

g(13)  =  l--^<0 

g{\4)  =  -!^ — 1<0 
U 

max 

0.1<;z,.<4.0  for  i  = 

where:  Uniax  =  0.125" 

Q  =  7Hz 
Syieid  =  78000  psi 
bi  =  0.5  hi 

hi  are  the  12  design  variables 
The  constraints  are  calculated  from  FEM  calculations 

The  constr  function  in  MATLAB  used  for  this  constrained  optimization  problem, 
uses  a  Sequential  Quadratic  Programming  (SQP)  method.  In  this  method,  a  Quadratic 
Programming  subproblem  is  solved  at  each  iteration,  and  an  estimate  of  the  Hessian  matrix 
of  the  Lagrangian  is  updated  at  each  iteration  using  the  BFGS  formula.  (Grace,  1992)  The 
termination  criterion  used  for  this  program  are  as  follows: 

Termination  criteria  for  h  =  IxlO'^ 

Termination  criteria  for  f  =  IxlO-^ 

Termination  criteria  for  g  =  IxlO^ 

Maximum  function  calls  =  1000 

The  program  was  run  for  6  different  initial  values  of  ho.  Table  5  shows  the  value 
of  the  objective  function  at  the  optimum  in  each  case.  Also  shown  are  the  values  of  the  14 


(17) 

(18) 

(19) 

(20) 


constraints  and  the  12  design  variables  corresponding  to  the  optimum.  The  last  row  in  the 
table  shows  the  number  of  function  evaluations  for  each  trial.  The  function  evaluation 
denotes  a  set  of  finite  element  analyses  of  the  frame  structure  for  both  static  and  eigenvalue 
analyses. 

The  results  show  that  the  design  was  improved  tremendously  over  the  initial 
design.  The  objective  function  for  the  4.0  initial  design  (worst  case)  was,  F(h®)  =  720.00 
in^.  Compared  with  the  corresponding  optimum  design  of  F(h*)  =  230.69  in^. 

From  the  results  in  Table  5,  it  can  be  seen  that  depending  on  the  starting  point,  the 
optimum  design  goes  towards  1  of  2  directions.  Initial  values  of  4.0,  3.9,  1.0,  and  0.1 
converge  towards  1  final  design.  While,  initial  values  of  3.0  and  2.0  converge  towards 
another  final  design.  The  different  designs  can  be  characterized  by  the  h(l)  and  h(12) 
values.  Both  designs  have  essentially  the  same  value  of  the  objective  function,  however 
the  designs  are  considerably  different.  It  appears  that  either  the  left  side  or  the  right  side  of 
the  frame,  but  not  both,  carries  the  bulk  of  the  load  depending  on  the  initial  design. 

An  additional  constraint  of  the  buckling  load  in  each  compressive  member  would  be 
useful  and  might  alleviate  the  problem  of  the  2  final  designs.  Also,  constraints  that  restrict 
the  relative  displacement  of  the  nodes  might  prove  useful.  Further  improvement  in  the 
convergence  of  the  optimal  designs  might  be  achieved  by  including  analytically  supplied 
gradient  information  to  the  optimizer,  such  as  partial  derivatives  of  the  objective  and 
constraint  functions. 
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Table  5.  Optimization  Results  for  6  Different  Initial  Design  Values 


Starting  Values  (ho)  for  all  12  Design  Variables  | 

4.0 

3.9 

3.0 

2.0 

1.0 

0.1 

230.6929 

230.7695 

230.6510 

230.9516 

230.9296 

233.8899 

g(i) 

-0.8938 

-0.8897 

-0.8671 

-0.8832 

-0.8886 

g(2) 

-0.8983 

-0.8973 

-0.8393 

-0.9035 

-0.9023 

g(3) 

-0.9062 

-0.8867 

-0.9058 

-0.9192 

g(4) 

-0.8771 

-0.8832 

-0.9085 

-0.8760 

-0.8938 

g(5) 

-0.9611 

-0.8775 

-0.8823 

-0.9622 

-0.9818 

g(6) 

-0.8801 

-0.8901 

-0.9572 

-0.9570 

-0.8946 

-0.9955 

g(7) 

-0.9817 

-0.9815 

-0.9701 

-0.9700 

-0.9834 

-0.9709 

g(8) 

-0.9664 

-0.9649 

-0.9844 

-0.9846 

-0.9723 

-0.9108 

g(9) 

-0.9110 

-0.9123 

-0.8792 

-0.9109 

-0.9820 

g(10) 

-0.9026 

-0.8945 

-0.9050 

-0.8896 

-0.9739 

g(n) 

-0.8369 

-0.8393 

-0.9488 

-0.9524 

-0.8381 

-0.8374 

g(12) 

-0.8696 

-0.8676 

-0.9404 

-0.9416 

-0.8669 

-0.8601 

g(13) 

-0.4135 

-0.4067 

-0.4060 

-0.4050 

-0.4055 

0.0005 

g(14) 

0.0025 

0.0001 

0.0037 

0.0001 

0.0008 

0.0000 

1.1282 

1.0585 

3.7739 

3.7482 

0.9787 

1.0441 

III 

1.0646 

1.0255 

2.4621 

2.4886 

1.0243 

1.0466 

3.1979 

3.2091 

1.1153 

1.0750 

3.1925 

Kifi 

2.0968 

2.1385 

1.5453 

1.5359 

2.1115 

2.7891 

ul 

2.1535 

2.1918 

1.4336 

1.4332 

2.2280 

0.3848 

h(6) 

1.4213 

1.4161 

2.2020 

2.1987 

1.4467 

0.3905 

h(7) 

2.8990 

2.9361 

2.3030 

2.3172 

2.9048 

3.4584 

h(8) 

2.3161 

2.3199 

2.8901 

2.9303 

2.3999 

2.3434 

h(9) 

1.5705 

1.5275 

2.1355 

2.1345 

1.5303 

0.2460 

h(10) 

1.1639 

1.0683 

3.1749 

3.2037 

1.0696 

0.1578 

h(ll) 

2.4631 

2.4830 

1.0405 

1.0222 

2.4995 

2.5613 

h(12) 

3.7649 

3.7474 

1.1347 

1.0588 

3.7331 

3.6312 

Active  Constraints 

14 

14 

14 

14 

14 

13,  14 

Function  Calls 

333 

340 

285 

395 

496 

857 
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2.  Solution  Using  Orthogonal  Array  Optimization  Technique 


The  same  problem  is  solved  again  using  the  technique  developed  in  this  research. 
The  optimization  problem  is  changed  slightly  for  implementation  of  this  technique.  Instead 
of  keeping  all  of  the  constraints  negative,  the  normalized  stress,  displacement,  and  lowest 
natural  frequency  are  maintained  between  0  and  1,  where  values  of  1  denote  active 
constraints. 

Therefore,  the  problem  statement  for  the  constrained  optimization  problem  is  given 
by  the  following  objective  function,  F(h),  and  constraints,  g(h): 


12 

minimize:  F(/i)  = 

i=\ 


(21) 


subject  to: 


g(0  = 


<1  for  i  =  1, 


Q. 

g(13)  =  — <1 
g(14)  =  -^<l 


12 


0.l<h.<4.0  for  i  =  1,...,12 


(22) 

(23) 

(24) 

(25) 


where.  ”  0.125 

£2  =  7Hz 
Syieid  =  78000  psi 
bi  =  0.5  hi 

hi  are  the  12  design  variables 

The  constraints  are  calculated  from  FEM  calculations 
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The  first  step  is  to  select  an  appropriate  orthogonal  array  to  deal  with  the  12  design 
variables  in  this  problem.  Therefore,  the  Li6(2i5)  Orthogonal  Array  (Taguchi,  1987)  as 
seen  in  Table  6  is  used. 

_ Table  6.  Li6(2^^)  Orthogonal  Array 


Ll6(2l5) 


Parameter  No. 

No.  of  Runs 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

1 

1 

1 

2 

2 

2 

2 

2 

2 

2 

3 

1 

1 

1 

2 

2 

2 

2 

1 

1 

1 

1 

2 

2 

2 

2 

4 

1 

1 

1 

2 

2 

2 

2 

2 

2 

2 

2 

1 

1 

1 

1 

5 

1 

2 

2 

1 

1 

2 

2 

1 

1 

2 

2 

1 

1 

2 

2 

6 

1 

2 

2 

1 

1 

2 

2 

2 

2 

1 

1 

2 

2 

1 

1 

7 

1 

2 

2 

2 

2 

1 

1 

1 

1 

2 

2 

2 

2 

1 

1 

8 

1 

2 

2 

2 

2 

1 

1 

2 

2 

1 

1 

1 

1 

2 

2 

9 

2 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

1 

2 

10 

2 

1 

2 

1 

2 

1 

2 

2 

1 

2 

1 

2 

1 

2 

1 

11 

2 

1 

2 

2 

1 

2 

1 

1 

2 

1 

2 

2 

1 

2 

1 

12 

2 

1 

2 

2 

1 

2 

1 

2 

1 

2 

1 

1 

2 

1 

2 

13 

2 

2 

1 

1 

2 

2 

1 

1 

2 

2 

1 

1 

2 

2 

1 

14 

2 

2 

1 

1 

2 

2 

1 

2 

1 

1 

2 

2 

1 

1 

2 

15 

2 

2 

1 

2 

1 

1 

2 

1 

2 

2 

1 

2 

1 

1 

2 

16 

2 

2 

1 

2 

1 

1 

2 

2 

1 

1 

2 

1 

2 

2 

1 

The  values  of  1  in  the  orthogonal  array  will  correspond  to  h  =  1.00",  and  the  values 
of  2  will  correspond  to  h  =  4.00".  The  FEM  programs  are  then  run  for  the  16  runs  in  the 
orthogonal  array  with  the  appropriate  values  of  h  for  each  design  variable.  Since  there  are 
only  12  design  variables,  the  last  3  columns  could  be  used  for  interactions,  but  will  not  be 
used  in  this  example.  The  results  from  the  running  of  the  orthogonal  array  are  shown  in 
Table  7. 
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s: 

<4-^ 
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VD 


<L> 

s: 


p 

CO 

Ji 

3 
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Weight 

45.00 

326.25 

326.25 

495.00 

382.50 

438.75 

438.75 

382.50 

382.50 

438.75 

438.75 

382.50 

382.50 

438.75 

438.75 

382.50 
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6.0302 
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Norm. 
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From  the  FEM  results  for  each  run  in  the  orthogonal  array,  the  averaged  Jacobian  is 
computed  (See  Equation  3)  which  contains  the  average  overall  sensitivities  of  each  of  the  12 
design  variables  to  the  14  constraints.  The  assumption  is  made  that  the  sensitivities  of  the 
normalized  stresses  in  each  member  due  to  other  members  is  small,  and  therefore  neglected, 
compared  to  the  sensitivities  of  the  normalized  stresses  due  to  itself.  Therefore  the  averaged 
Jacobian  matrix  is  simplified  as  shown  in  Table  8. 

The  absolute  maximum  value  in  each  column  is  then  used  to  determine  the  amount 
each  design  variable,  hi,  should  be  changed  because  the  absolute  maximum  value  is  the 


most  constraining.  It  is  denoted  by 


\^jj 


Equation  6  will  be  used  to  calculate  the 


value  Ahj,  where  Ahj  =  Axj . 

Table  9  shows  the  1-D  search  for  the  ho=[4  4444444444  4]  starting  point. 
The  highlighted  values  are  the  maximum  constraint  values  (1  is  the  target).  For  the  case 
when  Pi  =  2.2 ,  the  Maximum  Normalized  Displacement  exceeds  1,  therefore  the  distance  in 

the  search  direction  is  bracketed  and  simple  interpolation  is  used  to  get  a  value  close  to  1 . 

The  results  for  several  different  starting  points  are  shown  in  Table  10.  The  case 
when  ho  =  [111111111111]  is  initially  in  the  infeasible  domain,  with  a  resulting 
value  in  the  feasible  domain,  hi=[4  4444444444  4],  from  which  the  procedure  can  be 
repeated. 

Figures  17-20  show  the  displacements  of  the  frame  for  the  4  initial  designs.  The 
magnification  of  the  displacements  for  these  plots  is  50x.  Figures  21-23  show  the 
displacements  of  the  final  design  determined  by  using  the  orthogonal  array  for  3  of  the  initial 
conditions. 

Figures  24-26  show  the  first  3  mode  shapes  of  the  final  design  found  from  the  initial 
design  of  ho=[4  44444444444].  The  magnification  of  displacements  is  lOx. 
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where:  — -=  the  sensitivity  of  the  stress  in  element  1  to  height  1 

Ahi 


— the  sensitivity  of  the  maximum  displacement  to  height  i 
Ah- 


— —=  the  sensitivity  of  the  minimum  natural  frequency  to  height  i 
A/i, 


for  j  =  1  through  12 


33 


Table  9.  Optimization  Results  Using  Orthogonal  Array  with  ho=[4  4444444444  4] 


Function 

Evaluations 

Maximum 

Stress 

Constraint 

Minimum 

Normalized 

Natural 

Frequency 

Maximum  Weight 

Normalized 

Displacement 

16 

Array 

ho 

1 

1 

0.0677 

0.6804 

0.2898 

537.96 

P,.  =  1.2 

1 

0.0802 

0.7295 

0.3784 

471.93 

P.  =  1.4 

1 

0.0903 

0.7672 

0.4579 

430.50 

1 

-0.1348 

0.9144 

0.8748 

318.63 

P,.  =  2.2 

1 

-0.1573 

0.9795 

1.1195 

285.49 

p;.=2.i 

1 

0.1455 

0.9457 

0.9874 

301.80 

Figures  27-29  show  the  first  3  mode  shapes  of  the  final  design  found  from  the  initial 
design  of  ho=[3  3333333333  3].  The  magnification  of  displacements  is  lOx. 

Figures  30-32  show  the  first  3  mode  shapes  of  the  final  design  found  from  the  initial 
design  of  ho=[2  2222222222  2].  The  magnification  of  displacements  is  lOx. 

As  shown,  the  averaged  Jacobian  is  a  matrix  of  the  average  overall  sensitivities 
which  yields  a  search  direction  that,  when  coupled  with  a  1-D  search  technique,  can 
determine  a  design  which  lowers  the  objective  function  and  always  remains  in  the  feasible 
domain.  The  starting  point  value  is  updated  according  to 


X‘'  =  X^-'+aS‘‘ 


(26) 


where  Xq  is  the  new  vector  of  design  variables,  q  is  the  iteration  number,  S  is  a  search 
direction  vector  (here  calculated  by  the  averaged  Jacobian  matrix),  and  the  scalar  quantity  a 
is  the  distance  in  the  search  direction,  S. 
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Mesh  of  frame  for  12  elements 
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Figure  17.  Displacement  for  Initial  Design,  h=4.0.  Max.  Disp.  =  0.0195"  (Mag  =  50x) 

Mesh  of  frame  for  12  elements 

hs: 

3 
3 
3 
3 
3 
3 
3 
3 
3 
3 
3 
3 

0  10  20  30  40  50  60 

Xcoordiruite 


Figure  18.  Displacement  for  Initial  Design,  h=3.0.  Max.  Disp.  =  0.0595"  (Mag  =  50x) 
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Mesh  of  frame  for  1 2  elements 


X  coorcfinate 

Figure  19.  Displacement  for  Initial  Design,  h=2.0,  Max.  Disp.  =  0.2930"  (Mag  =  50x) 


Mesh  of  frame  for  12  elements 


0  .  10  20  30  40  50  60 
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Figure  20.  Displacement  for  Initial  Design,  h=1.0,  Max.  Disp.  =  4.6091"  (Mag  =  50x) 
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Mesh  of  frame  for  12  elemonls 
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Figure  21.  Displacement  for  Final  Design,  ho=4.0,  Max.  Disp.  =  0.1 185"  (Mag  =  50x) 

Mesh  of  frame  for  12  elements 
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Figure  22.  Displacement  for  Final  Design,  ho=3.0,  Max.  Disp.  =  0.1206"  (Mag  =  50x) 
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Mesh  of  frame  for  12  elements 
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Figure  23.  Displacement  for  Final  Design,  ho=2.0.  Max.  Disp.  =  0.1 185"  (Mag  =  50x) 
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Figure  24. 1st  Mode  Shape  with  ho=4.0  (Mag  =  lOx) 
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Second  Mode  Shape,  27.4  Hz 


First  Mode  Shape.  7.40  Hz 


Figure  25.  2nd  Mode  Shape  with  ho=4.0  (Mag  =  lOx) 

Third  Mode  Shape.  59.5  Hz 
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Figure  26.  3rd  Mode  Shape  with  ho=4.0  (Mag  =  lOx) 
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First  Mode  Shape.  7.56  Hz 


hs 

25361 

2.3985 

2.4965 

2.4206 

25427 

2.4956 

^361 8 

2.3633 

2.683 

2.6866 

25144 

Z6801 


Figure  27.  1st  Mode  Shape  with  ho=3.0  (Mag  =  lOx) 
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Figure  28. 2nd  Mode  Shape  with  ho=3.0  (Mag  =  lOx) 

Third  Mode  Shape.  5924  Hz 
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Figure  29.  3rd  Mode  Shape  with  ho=3.0  (Mag  =  lOx) 


-10  0  10  20  30  40  50  60 

Xeoordnals 


-10  0  10  20  30  40  50  60 

X  coordinate 


40 


First  Mods  Shape.  7.70  Hz 
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Figure  30. 1st  Mode  Shape  with  ho=2  (Mag  =  10) 
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Second  Mods  Shape,  2362  Hz 
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Figure  31.  2nd  Mode  Shape  with  ho=2  (Mag  =  10) 


h- 

2.4908 

2.6363 

26327 

2.6127 

2.4837 

26338 

^6751 

2.6625 

26353 

26315 

26138 

26384 


Figure  32.  3rd  Mode  Shape  with  ho=2,  (Mag  =  10) 


Third  Mode  Shape,  6068  Hz 
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Table  10.  Optimization  Results  (Orthogonal  Array)  for  4  Different  Initial  Designs 


Starting 

Values  (ho)  for  all  12  Design  Variables 

1 

i 

Q 

EBDuHl 

Objective  Function 

301.8008 

285.57 

284.45 

720.00 

g(l) 

0.1685 

0.0461 

g(2) 

0.1317 

0.0343 

g(3) 

0.0640 

0.0670 

0.0167 

g(4) 

ntiS 

0.0946 

0.0796 

0.0217 

g(5) 

-0.0148 

-0.0067 

-0.0075 

-0.0008 

g(6) 

-0.0173 

-0.0070 

-0.0071 

-0.0008 

g(7) 

-0.0182 

-0.0083 

-0.0084 

-0.0008 

g(8) 

-0.0168 

-0.0081 

-0.0086 

-0.0008 

g(9) 

-0.0764 

-0.0960 

-0.0212 

g(10) 

-0.0252 

-0.0488 

-0.0802 

-0.0162 

g(ll) 

-0.1001 

-0.1198 

-0.1317 

-0.0314 

g(12) 

-0.1455 

-0.1643 

-0.1762 

-0.0432 

g(13) 

0.9457 

0.9264 

0.9091 

0.5872 

g(14) 

0.9874 

1.005 

0.9878 

0.1629 

2.6261 

2.5361 

4.0 

2.2187 

2.3985 

2.6363 

4.0 

2.5088 

2.4965 

2.5327 

4.0 

h(4) 

2.2847 

2.4208 

2.6127 

4.0 

h(5) 

2.6458 

2.5427 

2.4837 

4.0 

h(6) 

2.5062 

2.4956 

2.5336 

4.0 

h(7) 

2.1100 

2.3618 

2.6751 

4.0 

h(8) 

2.1735 

2.3833 

2.6525 

4.0 

hC9) 

3.0613 

2.6830 

2.3353 

4.0 

E  1 

3.0719 

2.6866 

2.3315 

4.0 

2.5618 

2.5144 

2.5138 

4.0 

3.0526 

2.6801 

2.3384 

4.0 

Active  Constraints 

14 

14 

14 

N/A 

Function  Calls 

23 

22 

23 

17 

Final  Pj 

2.1 

1.15 

1.75 

0.9 

In  this  example,  a  simple  1-D  search  algorithm  is  used  which  extrapolates  out  in  the 
direction  chosen  until  a  constraint  is  violated.  Once  the  constraint  is  violated,  the  solution  is 
bounded  and  an  interpolation  procedure  similar  to  the  golden  section  method  may  be  used. 


3.  Comparison  of  the  two  Optimization  Techniques 


Table  11  shows  a  comparison  of  the  convergence  rates  between  the  SQP 
optimization  method  and  the  Orthogonal  Array  Technique.  It  is  important  to  notice  that  the 
SQP  routine  does  not  always  approach  the  optimum  from  the  feasible  domain.  However, 
the  Orthogonal  Array  Technique  always  approaches  the  optimum  solution  from  the  feasible 
domain.  This  is  an  extremely  important  feature  of  the  technique. 


Table  11. 

Convergence 

Comparison 

SQP 

Orthogonal  Array 

Function 

Objective 

Max 

Function 

Objective 

Max 

Evaluations 

Constraint 

Evaluations 

Constraint 

Value  <0 

Value  <1 

ho  =  {4.0} 

28 

444.596 

-0.5136 

23 

301.801 

.9874 

43 

-0.2585 

58 

274.131 

-0.554 

ho  =  {3.0} 

13 

405.000 

-0.2962 

22 

285.57 

1.005 

28 

300.724 

44 

283.249 

-0.1643 

ho  =  {2.0} 

13 

1.344  * 

23 

284.45 

0.9878 

28 

158 

1.945  * 

351 

q 

II 

o 

xa 

13 

45 

35.871  * 

17 

720 

0.5872 

496 

230.930 

.0008 

40  (17+23) 

301.801 

.9874 

*  indicates  infeasible  design 


For  the  ho  =  {2.0}  initial  condition,  the  SQP  converges  to  an  optimum  of  230.9  in 
351  function  evaluations  (with  infeasible  solutions  at  intermediate  steps),  and  the 
Orthogonal  Array  method  converges  to  284.5  in  only  23  function  evaluations.  The 
Orthogonal  Array  method  converges  to  within  23%  of  the  optimum  design  in  only  6.5%  of 
the  number  of  function  evaluations.  This  is  a  tremendous  computational  savings 
considering  the  FEM  analyses  required  at  each  function  evaluation.  Further  research  is 
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being  conducted  to  improve  the  convergence  closer  to  the  optimum,  but  studies  so  far  have 
been  inconclusive. 
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V.  CONCLUSION 


Improving  design  in  engineering  is  a  critical  process  in  order  for  businesses  to  stay 
competitive  and  profitable.  Therefore,  optimization  techniques  are  required  to  improve  the 
design  process  and  help  the  engineer  to  develop  better  products  in  the  future.  This  thesis 
has  described  a  new  design  optimization  technique  that  couples  the  orthogonal  array 
experimental  procedure  with  finite  element  analysis  in  an  effort  to  develop  an  economical 
and  efficient  way  of  optimizing  engineering  designs.  This  paper  has  also  shown  some 
examples  of  the  procedure  which  produced  very  good  results.  However,  the  proposed 
design  optimization  technique  needs  to  be  further  developed  and  refined  for  applications  to 
various  design  problems.  Results  also  show  that  the  orthogonal  array  technique  converged 
towards  the  optimum  at  a  faster  rate  than  the  SQP  method  for  the  frame  example. 
Additionally,  the  goal  of  reducing  the  number  of  FEM  calculations  was  achieved  as  well  as 
providing  the  engineer  with  insight  into  the  effects  of  changing  the  design  variables  on  the 
constraints.  This  sensitivity  information  may  prove  very  useful  when  developing  future 
designs  of  a  similar  nature. 
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VI.  RECOMMENDATIONS 


Further  coupling  of  the  orthogonal  array  technique  with  one-dimensional  search 
techniques  will  provide  for  a  more  robust  optimization  algorithm  and  should  be  the  topic  of 
future  research.  In  the  1-D  search  technique  a  possible  technique  might  be  to  vary  each 
value  of  the  target  value  individually,  in  an  effort  to  further  reduce  the  objective  function 
with  a  single  application  of  the  orthogonal  array. 

Also,  additional  studies  should  look  at  the  problem  of  refining  the  solution  further 
towards  the  actual  optimum.  One  possible  method  of  doing  this  might  be  to  incorporate  a 
second  orthogonal  array  around  a  refined  domain. 

Studies  should  be  conducted  on  additional  highly  non-linear  problems  to  further 
evaluate  the  performance  of  this  technique. 

The  cases  shown  in  this  thesis,  show  designs  that  mostly  start  from  the  feasible 
domain.  Additional  research  should  be  conducted  for  cases  of  an  initially  infeasible  design. 

Additionally,  utilizing  all  of  the  average  overall  sensitivities  and  not  just  the  diagonal 
terms  in  the  first  n  rows  of  the  averaged  Jacobian  should  be  examined  to  enhance  the 
procedure  developed  in  this  research. 

Finally,  research  conducted  bv  Yurkovich  (1994)  suggests  a  similar  approach  to  the 
optimization  of  structures  with  the  utilization  of  Taguchi  Methods.  However,  Yurkovich 
used  Taguchi  Methods  to  select  the  critical  design  variables  and  then  implimented  a 
conventional  optimization  technique. 

From  Yurkovich’s  research,  the  incorporation  of  interactions  into  the  procedure 
developed  in  this  thesis  might  prove  very  useful  in  improving  the  optimization  technique. 
Additionally,  orthogonal  arrays  can  be  further  utilized  to  reduce  the  number  of  design 
variables  and  interactions  that  must  be  analyzed  with  the  finite  element  method. 
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